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Abstract 

In this paper, we develop a multiscale model reduction technique that describes shale gas transport 
in fractured media. Due to the pore-scale heterogeneities and processes, we use upscaled models to 
describe the matrix. We follow our previous work [1], where we derived an upscaled model in the form 
of generalized nonlinear diffusion model to describe the effects of kerogen. To model the interaction 
between the matrix and the fractures, we use Generalized Multiscale Finite Element Method [SllIZ]- 
In this approach, the matrix and the fracture interaction is modeled via local multiscale basis functions. 
In [17], we developed the GMsFEM and applied for linear flows with horizontal or vertical fracture 
orientations on a Gartesian fine grid. In this paper, we consider arbitrary fracture orientations and use 
triangular fine grid and developed GMsFEM for nonlinear flows. Moreover, we develop online basis 
function strategies to adaptively improve the convergence. The number of multiscale basis functions in 
each coarse region represents the degrees of freedom needed to achieve a certain error threshold. Our 
approach is adaptive in a sense that the multiscale basis functions can be added in the regions of interest. 
Numerical results for two-dimensional problem are presented to demonstrate the efficiency of proposed 
approach. 


1 Introduction 

Shale gas transport is an active area of research due to a growing interest in producing natural gas from 
source rocks. The shale systems have added complexities due to the presence of organic matter, known as 
kerogen. The kerogen brings in new fluid storage and transport qualities to the shale. A number of authors, 
e.g., Loucks et al. (2009), Sondergeld et al. (2010), and Ambrose et al. (2012), |26l|3Tl|3], have previously 
discussed the physical properties of the kerogen using scanning electron microscopy (SEM) and showed the 
co-existence of nanoporous kerogen and microporous conventional inorganic rock materials. 

Gas transport in the kerogen typically develops at low Reynolds number and relatively high Knudsen 
number values. Under these conditions, it is expected that the transport is not driven by laminar (Darcy) 
flow dominantly but instead by the pore diffusion and other molecular transport mechanisms such as Knudsen 
diffusion and the adsorbed phase (or surface) diffusion. The latter introduces nonlinear processes at the pore 
scale that occur in heterogeneous pore geometry. Some types of upscaled models are needed to represent 
these complex processes for reservoir simulations. 

In large-scale simulations, the complex pore-scale transport needs to be coupled to the transport in frac¬ 
tures. This brings an additional difficulty in multiscale simulations. In particular, the multiscale simulations 
of the processes describing the interaction between the fracture and the matrix require reduced-order model 
approaches that work for problems without scale separation and high contrast. The objective of this paper 
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is to discuss the development of such approaches for describing the fracture and the matrix interaction by 
taking the upscaled matrix model following our previous work 

In our previous work [1] , we proposed a set of macroscopic models that take into account the nanoporous 
nature and nonlinear processes of the shale matrix. Our derivation uses multiple scale asymptotic analysis 
applied to mass balance equations, equation of state (for free gas) and isotherm of adsorption. The fine-scale 
microscopic description is largely based on the model formulated by Akkutlu and Fathi (2012), [2]. The 
macroscopic parameters that appear in the equations require solutions of cell problem defined in representa¬ 
tive volume elements (RVEs). These RVE problems take into account fine-scale variations and average their 
effects on macro scale. 

The multiscale approaches proposed in [T] are limited to representing the features that have scale sepa¬ 
ration. To represent the fracture network and the interaction between the fracture network and the matrix, 
we present a multiscale approach following the framework of Generalized Multiscale Finite Element Method 
(GMsFEM), [13] . The main idea of GMsFEM is to use multiscale basis functions to extract an essential 
information in each coarse grid (computational grid) and develop a reduced-order model. In [1^, we have 
developed the GMsFEM and applied for linear flows with horizontal or vertical fracture orientations on a 
Cartesian fine grid. In this paper, our contributions are: (1) the use of arbitrary fracture orientations and 
use triangular fine grids; (2) the development of GMsFEM for nonlinear flows; and (3) the development of 
online basis function strategies to adaptively improve the convergence. 

To represent the fractures on the fine grid, we use Discrete Fracture Model (DFM) |35|. The fine grid 
is constructed to resolve the fractures. For the coarse grid, we choose a rectangular grid. The GMsFEM 
framework uses these fine-scale models in computing the snapshot space and the offline space. The nonlinear 
models are handled with GMsFEM by locally updating multiscale basis functions. 

The study of flows in fractured media has a long history. Some modeling techniques on the fine grid 
include the Discrete Fracture Model (DFM), Embedded Fracture Model (EFM) [23 |25l [23], the single¬ 
permeability model, the multiple-permeability models (jSSlIllEHlESlEnilMlISQlfTOl), and hierarchical fracture 
models |23|. Though these approaches are designed for fine-scale simulations, a number of these approaches 
represent the fractures at a macroscopic level. For example, multiple-permeability models represent the 
network of connected fractures macroscopically by introducing several permeabilities in each block. The 
EFM ([23 [23 [23]) models the interaction of fractures with the fine-grid blocks separately for each block. 
The main idea of hierarchical fracture modeling presented in |23j is to homogenize small-length fractures 
(with the length smaller than the coarse block), while to represent the large-length fractures. Some of these 
approaches can be generalized by incorporating the interaction of fractures and permeability heterogeneities 
locally, which can lead to efficient upscaling techniques, mm- 

In recent papers EH, several multiscale approaches are proposed for representing the fracture effects. 
These approaches share common concepts with the methods that we discuss here in a sense that they add 
new degrees of freedom to represent the fractures on a coarse grid. The main difference is that our approaches 
use local spectral problems accompanied by adaptivity to detect the regions, where to add new basis functions. 
In this regard, the procedure of finding multiscale basis functions and the enrichment procedure is different 
from existing techniques. 

The proposed method constructs multiscale basis functions by appropriately selecting local snapshot space 
and the local spectral problems for the underlying nonlinear problem. The local spectral problems allow us 
to adaptively enrich in the regions with larger errors. In the paper, we discuss adaptivity issues and how to 
add multiscale basis functions in some selected regions. To reduce the computational cost associated with 
constructing the snapshot space, we follow [3 and use randomized boundary conditions. One of other novel 
components of the paper is the use of online basis functions (see [8] for online basis functions for steady state 
problems) for the time-dependent nonlinear problems. The online basis functions are constructed during the 
simulation using the residual and they can reduce the error significantly. These basis functions are used if 
the offline basis functions can not reduce the error below a desired threshold. 

We present numerical results for some representative examples. In these examples, we use nonlinear 
matrix and fracture models. Our numerical results show that the coarse-scale models with a fewer degrees 
of freedom can be used to get an accurate approximation of the fine-scale solution. In particular, only 10 % 
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degrees of the freedom are needed to obtain an accurate representation of the fine-scale solution. We also 
add a geomechanical contribution to the permeability term, where the permeability depends on the pressure. 
Furthermore, we demonstrate the use of online basis functions and how they can reduce the error. 

The paper is organized as follows. In the next section, we present a model problem. In Section 3, we 
discuss the fine-scale model. Section 4 is devoted to the development of GMsFEM, in particular, the offline 
spaces. In this section, we present numerical results for offline basis functions. In Section 5, we discuss 
randomized snapshot spaces and show that their use can give similar accuracy while for less computational 
cost. In Section 6 , we develop online basis functions and present numerical results. 


2 Model Problem 


In this paper, we will study nonlinear gas transport in fractured media motivated by several applications 
including shale gas. We are interested in the shale gas transport described in [2]. Similar equations arise in 
other models, where one considers a free gas in the tight reservoirs. We will consider general equations 

dc 

am(c)^ = div(&„(c,a;)Vc), (1) 

where c is the amount of free gas and a(c) and 6 (c) contain terms related to storage and adsorption coefficients. 
In [2], the authors consider the nonlinear terms have the forms 

dF OF Kj 

am{,c) = 0 ) 7 ^, bm{c,x) = (j)D + {1- + 4>-RTc, 

OC OC /i 

where 7 is a parameter, which is unity in kerogen and is equal to Vgrain,k/Vgrain in the inorganic material 
{Vgrain,k is grain volume and Vgrain is kerogen grain volume). Diffusivity D and porosity (j) are defined for 
the free fluid in the inorganic matrix and in the kerogen as follows 


D = 



in kerogen 
in inorganic matrix ’ 


in kerogen 
in inorganic matrix 


For the free gas we have ideal gas assumption. The Darcy law of free gas flow in inorganic matrix is 
used with permeability n and gas viscosity /r. For the sorbed gas we can use Langmuir or Henry‘s isotherms 
F = F{c). In |241135| . the authors discuss a general framework, where the equations also include nonlinear 
diffusivity due to adsorbed gas in a shale formation. In [22] . the nonlinear terms appear due to barotropic 
effects. 

The nonlinear flows also contain components that are due to diffusion in the fractures. One needs 
additional equations for modeling fractures. The fractures have high conductivity. We will use a general 
equation of the form 

dc 

afic)—= diY{bf{c,x)Vc) (2) 

to describe the flow within fractures. In [5], the authors use 

af{c) = (j)f, bf{c,x) = —RTc, 

where 0/ and k/ are the fracture porosity and permeability. These problems are solved on a fine grid using 
DFM as will be described in Section |3l 

In many shale gas examples, the matrix heterogeneities can be upscaled and the resulting upscaled equa¬ 
tion has the form 0 . However, the interaction between matrix and fractures require some type of multiscale 
modeling approach, where the effects of the fractures need to be captured more accurately. Approaches, such 
as multicontinuum [23] are often used, but these approaches use idealized assumptions on fracture distribu¬ 
tions. In this paper, we will use multiscale basis functions to represent fracture effects. In our previous work 
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m, we have considered similar approaches for single-phase flow when fractures (which could be horizontal 
or vertical) are aligned with Cartesian grid. In this paper, we consider arbitrary fracture distribution in the 
context of nonlinear flow equations. 

The overall model equations will be solved on a coarse grid. Next, we introduce the concepts of fine and 
coarse grids. Let be a usual conforming partition of the computational domain into finite elements 
(triangles, quadrilaterals, tetrahedra, etc.). We refer to this partition as the coarse grid and assume that 
each coarse element is partitioned into a connected union of fine grid blocks. The fine grid partition will be 
denoted by T^, and is by definition a refinement of the coarse grid ■ We use {xi}^i (where N denotes 
the number of coarse nodes) to denote the vertices of the coarse mesh and define the neighborhood of 
the node Xi by 

^ S (3) 

See Figure for an illustration of neighborhoods and elements subordinated to the coarse discretization. 


(Coarse Grid) 



Figure 1: Illustration of a coarse neighborhood and coarse element 

We emphasize the use of to denote a coarse neighborhood, and K to denote a coarse element throughout 
the paper. 

3 Fine-scale discretization 

To discretize the system on fine grid, we will use finite element method and use DFM for fractures. To 
solve Problem Q using finite element method (FEM), we need a fine grid discretization to capture the 
fractures. These computations can be expensive. Here, we apply the discrete fracture network (DFM) model 
for modeling flows in fractures [29]. 

In the discrete-fracture model, the aperture of the fracture appears as a factor in front of the one 
dimensional integral for the consistency of the integral form. This is the main idea of the discrete-fracture 
model, which can be applied in any complex configuration for fractured porous media. 

To demonstrate it, we consider the two-dimensional problem of Equation ([^. We simplify the fractures as 
the lines with small aperture. Thus, one-dimensional element is needed to describe fractures in the discrete- 
fracture model. The system of equations Q will be discretized in a two-dimensional form for the matrix 
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and in one-dimensional form for the fractures. The whole domain can be represented by 




( 4 ) 


where m and / represent the matrix and the fracture of the permeability field k, respectively. Here, di is the 
aperture of the i th fracture and i is the index of the fractures. Note that is a two-dimensional domain 
and idf i is a one-dimensional domain (Figure]^. Then Equations Q and ([^ and can be written as follows 
(for any test function v): 


m{ 




v) -I- a(c, v) 




amic)—vdx + di'Y^ J a/(c)—uda;-|- 

bm(c,x)'Vc-'Vvdx + di2^ / 6m(c, a:)Vc • Vu dx = 0. 

Jn f . 


(5) 



Figure 2: Fine grid with fractures 


To solve ([^, we will first linearize the system. We will use the following linearization 

„n+l _ f „n+l _ n _ r „n+l _ n 

m{ - ,v) + ,v) = / amic^) - vdx + di^^ / af(c^) - vdx+ 

. Jcifi 


+ [ bjn{c^,x)\7c'^'^^ ■ \7v dx + di^^ ( 6/(c", a;)Vc"'''"^ • Vi; dx = 0. 

J Otti „■ J i 


( 6 ) 


The standard fully-implicit finite difference scheme is used for the approximation with time step size r and 
superscripts n, n-|-1 denote previous and current time levels. This is a first-order in time and unconditionally 
stable linearization. 

For standard Galerkin finite element method, we write the solution as c = where (pi are the 

standard linear element basis functions defined on and Nf denotes the number of the nodes on the fine 
grid. The equation § can be presented in matrix form: 


fXl + l 




— c 


-A 


n n+1 


= 0 , 


(7) 
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where M is the mass matrix given by 


M” = [wy] = / am{c^)4'i(l)j dx + di'^ / af{c'^)(j)r(j)jdx, 

J Orn 2 ^ f 

and A is the stiffness matrix given by 

A'^ = [aij] = / bm{c^, x)^/4>i ■ \/(pj dx + di^^ / b f , x)^/4>i ■(pj dx. 

J ^ J Vlf 

Hence at each time step we have the following linear problem 

g«^n+l ^ 

where Q" = (M” + tA^). This fine scale discretization yields large matrices of size Nf x Nf. 


( 8 ) 


4 Coarse-grid discretization using GMsFEM. Offline spaces. 

We use multiscale basis functions to represent the solution space. We will consider the continuous Galerkin 
(CG) formulation and signify uji as the support of basis functions. We denote the basis functions by ip '^', 
which is supported in uJi, and the index k represents the numbering of these basis functions. In turn, the 
GG solution will be sought as 

Cms{x,t) = ^cl{t)xPl'{x). 

i^k 

Once the basis functions are identified, the GG global coupling is given through the variational form 

Oc 

+ a{cms,v) = 0, for all u e Goff, (9) 

ot 

where Vbff is used to denote the space spanned by those basis functions and 


m{c,v) = / amCvdx + di'S^ / afcvdx, 
i(c,v)= / bmS c ■ V V dx + di'S^ / bfWc-S/vdx. 

J ■ J Cl f i 


Let V be the conforming finite element space with respect to the fine-scale partition T^. We assume 
c G G is the fine-scale solution satisfying 


Oc 

m( — ,v) + a(c,v) = 0, vGV. 
ot 


( 10 ) 


Next, we describe GMsFEM. GMsFEM consists of offline and online stage. In the offline stage we 
construct multiscale basis functions and after that in the online stage, we solve our problem for any input 
parameters, such as right hand sides or boundary conditions. 

Offline computations: 


Step 1. Coarse grid generation. 

Step 2. Construction of the snapshot space that will be used to compute an offline space. 


Step 3. Construction of a “small” dimensional offline space by performing dimension reduction in the 
space of local snapshots. 
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Given the computational domain, a coarse grid can be constructed and local problems are solved on 
coarse neighborhoods to obtain the snapshot spaces. Then, smaller dimensional offline spaces are obtained 
from the snapshot spaces by dimension reduction via some spectral problems [m [121 mi mi m. After that 
we can solve our problem in the constructed offline space. Moreover, we will construct online basis functions 
that are problem dependent and are computed locally based on some local residuals Pg. 

We now present the construction of the offline basis functions and the corresponding spectral problems 
for obtaining a space reduction. In the offline computation, we first construct a snapshot space I^nap- The 
snapshot space can be the space of all fine-scale basis functions or the solutions of some local problems with 
various choices of boundary conditions. For example, we can use the following K-harmonic extensions to form 
a snapshot space. For each fine-grid function, 5^{x), which is defined by 5^{x) = Sj^k, G Jh{uJi), where 
Jh{uJi) denotes the fine-grid boundary node on duJi. For simplicity, we omit the index i. Given a fine-scale 
piecewise linear function defined on duj (here cu is a generic coarse element), we define by following 

variational problem 

, v) = f -Vv dxdj'^^ f b fV ■ V v dx = Q in w, (11) 

Juj^ j Jujfj 

and = Sj{x) on 9a;, co = uJm ©j djojfj. 

For brevity of notation, we now omit the superscript oj, yet it is assumed throughout this section that the 
offline space computations are localized to respective coarse subdomains. Let li be the number of functions 
in the snapshot space in the region to, and 

Fsnap = span{'0™‘*''’ : 1 < j < k}, 


for each coarse subdomain ui. 

Denote 

In order to construct the offline space we perform a dimension reduction of the snapshot space using 
an auxiliary spectral decomposition. The analysis in m motivates the following eigenvalue problem in the 
space of snapshots: 

^ (^ 2 ) 

where 

= [a±] = f dx + d, ^ dx = 

JuJm j 

= [c® ] = [ dx+d,Y,[ dx = 

Ju:m j 

where A and S denote analogous fine scale matrices as defined by 

Aij = / bm(c^ ,x)V(l)i(j)j dx + di'S^ / bf{c^,x)V(pi •V(pj dx, 

JOrr, , JDf_, 

Sij — / 6771(0 dx © di ^ ( j bj^i^c dx, 

7 JOf,. 

where (j)i is the hne-scale basis function. To generate the offline space, we then choose the smallest 
eigenvalues from Eq. ( |12[ ) and form the corresponding eigenvectors in the space of snapshots by setting 
(fop k = 1 ,..., M^g), where are the coordinates of the vector 'I'))®. 

Next, we create an appropriate solution space and variational formulation that for a continuous Galerkin 
approximation. We begin with an initial coarse space = span{xi}N;^. Recall that N denotes the 
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number of coarse neighborhoods. Here, Xi ^re the standard multiscale partition of unity functions defined 

by 


• Vu dx + dj / 6/Vxi • Vu dx = 0 K G oj 

Xi = gi on dK, 


(13) 


for all K G u>, where gi is a continuous function on dK and is linear on each edge of dK. 

We then multiply the partition of unity functions by the eigenfunctions in the offline space 
construct the resulting basis functions 


to 


= Xi'^T'°^ for 1 < j < TV and 1 < A: < 


(14) 


where denotes the number of offline eigenvectors that are chosen for each coarse node i. We note that 
the construction in Eq. (14) yields continuous basis functions due to the multiplication of offline eigenvectors 


with the initial (continuous) partition of unity. Next, we define the continuous Galerkin spectral multiscale 
space as 

VoS = span{'!/>i^fe : 1 < i < TV and 1 < A: < M^^}. (15) 

Using a single index notation, we may write VoS = span{'0i}^i, where TVc = denotes the total 

number of basis functions in the space VoS- We also construct an operator matrix 

Rq = 

where ipi are used to denote the nodal values of each basis function defined on the fine grid. 

We seek c„is(a^) = Ciipi{x) G UoS such that 




’j) + a(cnis, v) = 0 for all v G Uoff- 


We note that variational form in (161 yields the following linear algebraic system 




0 ^0 5 


(16) 


(17) 


where cq denotes the nodal values of the discrete CG solution, and Qq = RoQ’^Rq and Mq = RqM^. We 
also note that the operator matrix may be analogously used in order to project coarse scale solutions onto 
the fine grid In our simulations presented next, we do not update basis functions. We 

discuss basis function update in Section]^ 


4.1 Numerical result 

We present numerical results for the coarse-scale solution using offline basis functions. The basis functions 
of the offline space are constructed following the procedure described above. Note that, the basis functions 
are constructed only once at initial time and used for generating the stiffness matrix and the right hand side. 

We consider the solution of problem with constant and nonlinear matrix-fracture coefficients in ([^. As 
constant coefficients (see previous section) representing matrix and fracture properties, we use following 

Om = 0.8, 6m = 1-3 • 10“^ and a/= 0.001, 6/= 1.0. 


For nonlinear matrix-fracture coefficients, we use 


dF 

am{c) = 


d F Kj 

bmiCjx) = (l)D + {l- (j))Ds— + (j)-RTc, 


(18) 


af{c) = cj)f, bf{c,x) = —RTc, 


and 


( 19 ) 





where Dk = 10 Di = 10 ®[m^/s], (j) = 0.04, T = 413[i^], ^ = 2 • 10 ^[kg/{ms)] and for fractures 

kf = 10-i2[m2], (jjf = 0.001. 


As for permeability k in (18), we use constant k = kq and stress-dependent model k = Km (see mm) 
with 


M'^ 




Ko 1 - 


Pc - ap 


Pi 


where Kq = 10 p = RTc, Pc 

gas, we use Langmuir model 


10®[Po], Pi = 1.8 • 10®[Po], a = 0.5 and M = 0.5. For the sorbed 


F{c) = c^s 


s 

(1+SC)2’ 


where s = 0.26 • 10 ® and c^s = 0.25 • 10 ^[mol/m^]. 



Figure 3: Coarse and fine grids. Coarse grid contains 50 cells, 85 facets and 36 vertices. Fine grid contains 
7580 cells, 11470 facets and 3891 vertices. 




Figure 4: Coarse and fine grids. Coarse grid contains 200 cells, 320 facets and 121 vertices. Fine grid contains 
13036 cells, 19694 facets and 6659 vertices. 

The equation is solved with Dirichlet boundary condition c{x, t) = 5000 on the left boundary and Neu¬ 
mann boundary conditions = 0 on other boundaries. The domain fl has a length of 60 meters in both 
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directions. We calculate concentration for t^ax = 5 years with the time step r = 10 days. As for initial 
condition, we use c{x,t = 0) = 10QQ0[mol/m^]. For the numerical solution, we construct structured two 
coarse grids with 36 nodes (Figureand with 121 nodes (Figure]^. As for fine grids, we use unstructured 
grids, which resolves the existing fractures. 


Fine 


Coarse 



Figure 5: Solution with constant matrix-fracture coefficients on coarse (top) and on fine (bottom) grids for 
t=l, 3 and 5 year (from top to bottom) 

In Figure we show the pressure distribution for three concrete time level t = 1, 3 and 5 years. For 
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Figure 6: Coarse-scale (top) and fine-scale (bottom) solutions for t=5 year for the case of nonlinear perme¬ 
ability with K = kq 



Figure 7: Coarse-scale (top) and fine-scale (bottom) solutions for t=5 year for the case of nonlinear coefficients 
with K = Km 


the pressure and concentration, we have the following relationship: p = RTc. Pressure distribution for 
nonlinear matrix-fracture coefficients in Q is presented in Figures |^-[^ for last time level. In these figures, 
we show fine-scale (reference) and coarse-scale (multiscale) solutions. The coarse-scale solution is obtained 
in an offline space of dimension 288 (using Mof / = 8 multiscale basis functions per coarse neighborhood) and 
the fine-scale solution is obtained in a space of dimension 3891. Compared to the fine-scale solution on the 
left with the coarse-scale solution on the right of the figures, we observe that the GMsFEM can approximate 
the fine-scale solution accurately. 

To compare the results, we use relative weighted errors 

11 ^ 11 * II ^ms 11 * /11 II*; 
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Figure 8: Multiscale basis functions corresponding to the first 3 smallest eigenvalues in the case with constant 
fracture-matrix properties after multiplication to partition of unity functions, i = 25 and 

A: = 0,1, 2 (from left to right) 


In Table b we present relative errors (in percentage) for last time level for constant fracture and matrix 
properties in using coarse grids with 36 and 121 nodes. For the coarse-scale approximation, we vary the 
dimension of the spaces by selecting a certain number of offline basis functions {Moff) corresponding to the 
smallest eigenvalues. In the Table we recall that 14// denotes the offline space, dim{Voff ) is the offline 
space dimension, M^ff is the number of the multiscale basis functions per coarse neighborhood (we use a 
similar number of M^ff for each w^), Cms and c/j are the multiscale and reference solutions, respectively. 

Figurej^presents the multiscale basis functions corresponding to the first 3 smallest eigenvalues in the case 
with constant fracture-matrix properties in Q. These offline basis functions are multiplied by partition of 
unity functions. When we use Moff = 8 and the case with 36 coarse nodes, the relative and weighted 
errors are 0.3 % and 0.7%, respectively for hnal time level. The dimension of the corresponding offline space 
is 288 and for reference solution is 3891. For coarse grid with 121 nodes, the relative errors are slightly 
smaller 0.1% and 0.2% for and weighted errors, respectively. The dimension of the corresponding 
offline space is 968 and for reference solution is 6659. The relative and errors at different time instants 
for the cases with 36 and 121 coarse grids are presented in Figures and IT As we observe if we take 4 or 
more basis functions per coarse node, the relative errors remain small. 


Moff 

dim(t4//) 

^min 

T ^ 

HL 

1 

36 

9.0 10-y 

24.484 

84.383 

2 

72 

4.5 10-® 

12.229 

33.923 

4 

144 

1.1 10-'^ 

1.068 

2.162 

8 

288 

2.2 10-6 

0.303 

0.737 

12 

432 

0.19 

0.083 

0.258 


Moff 

dim(t4//) 

^min 



1 

121 

2.5 10-6 

17.136 

68.989 

2 

242 

9.6 10-8 

3.975 

36.337 

4 

484 

1.6 10-'^ 

0.651 

3.595 

8 

968 

0.37 

0.110 

0.246 

12 

1452 

1.23 

0.060 

0.108 


Table 1: Numerical results (relative errors (%) for the final time level). Left: for the case with 36 coarse 
nodes. Right: for the case with 121 coarse nodes. 


We present relative weighted errors in Tables and for different number of eigenvectors Mof / for the 
case with nonlinear matrix-fracture coefficients in §. We consider a case with 36 coarse nodes. When we 
use Moff = 8 and the case with k = kq, the relative and errors are 0.2 % and 0.7%, respectively. 
The dimension of the corresponding offline space is 288 and for reference solution is 3891. For the case with 
K = Km in (18), we have 0.4% and 1.0% of relative and errors, respectively. The dimension of coarse 
spaces for the corresponding number of eigenvectors are 72, 144, 288, 432 and 576 for Moff = 2,4,8 and 12. 
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Figure 9: Relative and H\ weighted errors (%) for coarse grid in Figure with 36 nodes. Constant 
matrix-fracture properties. 
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Figure 10: Relative and weighted errors (%) for coarse grid in Figure]^ with 121 nodes. Constant 
matrix-fracture properties. 



We observe that as the dimension of the coarse space (the number of selected eigenvectors Moff) increases, 
the respective relative errors decrease. Also we have similar error behaviour as for case with constant matrix- 
fracture coefficients. Moreover, we see that the decrease in the relative error is fast initially and one can 
obtain small errors using only a few basis functions. 


^off 

dim{Voff) 

^min 

r ‘2 




dim(Vo//) 

^min 

r‘2 

hT 

1 

36 

4.8 10"® 

21.717 

87.897 


1 

121 

2.5 10"® 

14.333 

64.197 

2 

72 

2.4 10-® 

10.772 

38.774 


2 

242 

9.6 10-® 

3.673 

30.510 

4 

144 

6.0 10-® 

0.933 

1.947 


4 

484 

1.6 10-^ 

0.646 

3.272 

8 

288 

1.1 10"® 

0.270 

0.737 


8 

968 

0.37 

0.110 

0.251 

12 

432 

0.19 

0.123 

0.323 


12 

1452 

1.23 

0.063 

0.159 


Table 2: Numerical results (relative weighted errors (%) for final time level) for case with k 
Left: the case with 36 coarse nodes. Right: the case with 121 coarse nodes. 


kq in (181. 
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^Off 

dim(Vo//) 

^min 

r 2 


1 

36 

2.0 10-y 

24.484 

92.039 

2 

72 

1.0 10"® 

10.785 

35.874 

4 

144 

2.6 10"® 

1.247 

2.423 

8 

288 

5.1 10-^ 

0.432 

1.069 

12 

432 

0.19 

0.234 

0.712 


Moff 

dim(Vo//) 

^min 

r‘2 


1 

121 

2.5 10"® 

16.318 

60.267 

2 

242 

9.6 10-® 

3.715 

22.872 

4 

484 

1.6 10-^ 

0.645 

3.000 

8 

968 

0.37 

0.134 

0.386 

12 

1452 

1.23 

0.096 

0.282 


Table 3: Numerical results (relative weighted errors (%) for the final time level) for case with k = Km in 
(18). Left: the case with 36 coarse nodes. Right: the case with 121 coarse nodes. 


Remark 4.1. In our numerical simulations, we do not use empricial interpolation procedures for approx¬ 
imating the nonlinear functionals a.(c,-) and b.(c,-) (see for more details). In the approaches of 
empirical interpolation concepts are used to evaluate the nonlinear functions by dividing the computation 
of the nonlinear function into coarse regions, evaluating the contributions of nonlinear functions in each 
coarse region taking advantage of a reduced-order representation of the solution. By using these approaches, 
we can reduce the computational cost associated with evaluating the nonlinear functions and consequently 
making the computational cost to be independent of the fine grid. 


5 Randomized oversampling GMsFEM 


Next, we present numerical results for the oversampling and the randomized snapshots that can substantially 
save the computational cost for snapshot calculations. In this algorithm, instead of solving local harmonic 
problems 0 for each fine grid node on the boundary, we solve a small number of harmonic extension local 
problems with random boundary conditions [S]. More precisely, we let 

/cjj.rsnap ‘T + 

iPj =rj, X £ OLof , 


where rj are independent identical distributed standard Gaussian random vectors on the fine grid nodes of 
the boundary. When we use randomized snapshots, we only generate a fraction of the snapshot vectors by 
using random boundary conditions. 

For snapshot space calculations, we use the extended coarse grid neighborhood for m = 1,2,..., by 
w)*’ = oji -\- m, where m is width of the fine-grid layer. Here, for example, ujf' = oji -\- 1 means the coarse 
grid neighborhood plus all 1 layer of adjacent fine grid of Wi, and so on (see Figure 11 for illustration). 
Calculations in the oversampled neighborhood domain ujf reduces the effects due to the artificial oscillation 
in random boundary conditions. 


5.1 Numerical results 

The simulation results are presented in Tables i-0 for 36 node coarse grid case. We use constant matrix- 
fracture properties, see (|^. We present the results for the randomized snapshot case for last time level. In 
our simulations, we set the oversampling size m = 0, 2,4,6 for ojf = oJi -\-m and use different numbers of 
multiscale basis functions M^ff = 2,4,8 and 12. 

In Table we investigate the effects of the oversampling wj*" = uJi-\-m, as we increase the number of fine 
grid extensions m = 0, 2,4 and 6. We see that the oversampling helps to improve the results initially, but 
the improvements slow and larger oversampling domains do not give significant improvement in the solution 
accuracy. When we use a snapshot ratio of 25.6 % (between the standard number of snapshots and the 
randomized algorithm for wj*" = uji -\-4), the relative and H) weighted errors are 0.2 % and 0.8% for full 
snapshots and 0.2 % and 0.9% for randomized snapshots. We observe that the randomized algorithm can 
give similar errors as a full snapshots. 
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Figure 11: Neighborhood domain with oversampling = Wi + m, m = 1,2, 4, 6) for the coarse grid with 
36 nodes 

Tableshows relative and errors for different number of randomized snapshots Mi. The oversam¬ 
pled region uj^ = -I- 4 is chosen, that is, the oversampled region contains an extra 4 fine-grid cell layers 

around Wi. Our numerical results show that one can achieve a similar accuracy when using a fraction of 
snapshots with randomized algorithms and thus, it can provide a substantial CPU savings. 


6 Residual based adaptive online GMsFEM 


In this section, we consider the construction of the online basis functions that are used in some regions 
adaptively to reduce the error significantly. We follow earlier works mi, which were done for linear time- 
independent problems. The online basis functions are constructed based on a residual and take into account 
distant effects. The construction of online basis functions is motivated by the analysis. Using the offline 
computation, we construct multiscale basis functions that can be used for any input parameters to solve 
the problem on the coarse grid. The fast convergence due to adding online basis functions depends on the 
offline space. It is important that the offline space contains some essential features of the solution space. In 
our numerical simulations, we demonstrate that with a sufficient number of offline basis functions, we can 
achieve a rapid convergence for the proposed online procedure. 

First, we derive the error indicator for the error (c" — for time-dependent problem (21) in the energy 


norm. Furthermore, we use the error indicator to develop an enrichment algorithm. The error indicator gives 
an estimate of the local error on the coarse grid region Wi and we can then add basis functions to improve 
the solution. 

We assume, as before, V is the fine-scale finite element space. To find the fine-scale solution G V, 
we solve (as before) 

n+l _ n 

,u)-f a(c”+\u) = (/,u). 


TO 


MvGV 


and for multiscale solution S Voff we have 


n+l _ „n 

mC-^ -^,^) + a(cS,r,u) = (/,u), VueUoff. 


( 20 ) 


( 21 ) 


We define a linear functional r"(u) for n-th time level by 


r"(u) =t(/,u) -TO(c”t^ -™(cmt\t')- 
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Moff 

full snapshots 

Ll Hi 

randomized snapshots 

Ll Hi 

without oversampling, uJi 


100 % 

39.7 % 

2 

12.229 

33.923 

8.303 

33.237 

4 

1.068 

2.162 

1.704 

4.730 

8 

0.303 

0.737 

1.005 

2.962 

12 

0.083 

0.258 

0.557 

1.643 


with oversampling, = 

- U)i 2 


100 % 

28.3 % 

2 

12.247 

33.943 

8.921 

33.399 

4 

1.073 

4.237 

0.972 

3.750 

8 

0.261 

0.744 

0.354 

1.003 

12 

0.114 

0.329 

0.219 

0.704 


with oversampling, = 

- LOi + 4 


100 % 

25.6 % 

2 

12.216 

33.657 

9.334 

28.213 

4 

1.015 

4.576 

0.626 

2.561 

8 

0.262 

0.841 

0.264 

0.949 

12 

0.114 

0.349 

0.153 

0.441 


with oversampling, = 

= coi 6 


100 % 

22.5 % 

2 

12.746 

35.899 

9.455 

27.922 

4 

1.013 

5.014 

0.603 

2.377 

8 

0.251 

0.820 

0.277 

0.875 

12 

0.124 

0.369 

0.120 

0.421 


Table 4: Randomized oversampling for GMsFEM with number of snapshots Mi = 24 (constant matrix- 
fracture properties) in every = cui + n, n = 0,2,4,6 for coarse mesh with 36 nodes (relative errors (%) 
for final time level) 


LLoff 

12.8 % {M, = 12) 
i/i 

17.0 % (M, = 16) 

02 i/1 

21.3 % {Mi = 20) 
02 

25.6 % (M, = 24) 

02 071 

29.8 % (M, = 28) 
02 ijl 

2 

8.228 

24.878 

9.449 

28.895 

7.346 

22.774 

9.334 

28.213 

9.335 

26.973 

4 

1.908 

4.208 

1.381 

3.581 

0.779 

2.692 

0.626 

2.561 

0.843 

4.439 

8 

0.861 

1.777 

0.589 

1.563 

0.292 

1.189 

0.264 

0.949 

0.245 

0.894 

12 

- 

- 

0.313 

0.781 

0.217 

0.581 

0.153 

0.441 

0.110 

0.393 


Table 5: Randomized oversampling for GMsFEM with different number of snapshots Mi = 12,16,20,24,28 
in every ojf = cvi + 4 (constant matrix-fracture properties) for coarse mesh with 36 nodes (relative errors 
(%) for final time level) 


’ dx 


Let Wi be a coarse region and V) = Hg^uji) then 

= r f fv- f am (c”t^ “ c”s) v dx - t f bm Vc”t^ ' Vv < 

J OJi J J 

“ H If - c”s) vdx-rdjY^ f bf Vc”t^ • Vu dx, 

J j 

where uji = uj'^ djUj{’^ and m and / represent the matrix and the fracture. 
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The solution at (n + 1) time level is the solution of the elliptic problem of the form 


ar(Ci ,v)=T{f,v)+m{c^s,v). 


(22) 


We use following notation 

ar{u, v) = m(u, v) + Ta{u, v). 

Error estimators for the spatial discretization error take into account the dependence of the elliptic problem 


(22) on the time step parameter t and we will use the t- weighted Hi norm 




where 

\\v\\l = ar{v,v), \\v\\l=a{v,v), \\v\\l^ = m{v,v). 

We define the projection 11 : y Vojf by 

N 

where Pi : V ^ spanji/)^’’"-^-^} be the projection defined by 


h 

PiV = ^ ■ 

k=l 

The projection Pi is the first li terms of spectral expansion in terms of eigenfunctions of following problem 


(vI/o«,^) = A°«r 6|Vx 


P 'i'°^vdx. 


(23) 


Then 


and 


■ [ b\\/Xz\‘^ {v - P^vf dx < - PiV,v- Piv). 

J (jJi ^h+1 


(?; - PiV, V - Piv) < (?;, v). 


We note that this spectral problem is different from the original one formulated in (12); however, it involves 
similar terms, such as energy norms and norms. 


Let e" = c" — c”s is error for n-th time level and using (20) and (21), we have 

m(e”+i - e", v) + ra(e”+\ i;) = (r”, v), 


(24) 


where the right hand side can be written as follows 

r"(r;) = T{f,v) - m(c"t^ - - ra{c^^,v) 

< r{f, r; - n?;) + r(/. Hi;) - - c”„ Hz;) - ro(c”t\ Hz;) 

- - c”s> - nz;) - Ta{c^\v - Hz;) 

= r{f,v- Hz;) - m{c^^ - c”„z; - Hz;) - Ta(c”t^^' “ Hz;) 

AT N 

= H “ PiV)\\r- 

2=1 2=1 
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We have 


\\x^iv - P^v)\\l = \\x^iv - PivWm + rWxiiv - Piv)\\l 

= WXtiv - Piv)\\l, + tC f b\VXtl"^ (v - P^vf dx+ tC f b(v - P,v)\‘^ dx 
J OJi J OJi 

< Ca‘^'{v — PiV,v — Piv) + tC f b |Vxip {v — Piv)^ dx 

J OJi 


(25) 


< [C + {v - PiV, V - Piv) < ( C + 


c 

^k+i 


aT'iv,v). 


Therefore 


iV JV ^ „ V i/Z 

(e"+i - e\v)+Taie-+\v) < ^ ||rriU||x.(^ - P^v)\\r < E H^riU + —j a-'(v,v)P^ 

N \ 1/2 AT ^ X 1/2 / Af 

Eii’-riK i: s (c + —) E 

i=l / i=l ^ min/ \i=l 


< c + 


1/2 / N 


. 1/2 

z”ll*) ariv, 


where A^in = min* . 
Finally, we take v = 


n+ 1 112 


i< c + 


c 


A„ 


1/2 / N 


El 


< ic + 


c 




1/2 / AT 


1/2 


1/2 


El 


^n\\2 


o'^+l I 


n+ 1 1 


, + m(e",e’^+i) 


|e"ILI|e"+i| 


Then 


|e-+i||, < [C + 


C 


A„ 


1/2 / N 


1/2 


El 


i+ii? 


This inequality residuals give a computable indicator of the error e"''"! = in the r-weighted Hi 

norm. 


Remark 6.1. ITe note that the analysis suggests the use of (23) as a local eigenvalue problem. This 
eigenvalue problem is “slightly” different from (12) that we have used earlier. Our numerical simulations 
show that the use of (23) improves the convergence of the offline or online procedures slightly in our numerical 


examples. We will use the spectral problems based on in our numerical simulations as it is independent 
of time stepping. 


Next, we consider online basis construction. We use the index m > 1 to represent the enrichment level. 
At the enrichment level m, we use Vff. to denote the corresponding space that can contains both offline and 
online basis functions. We will consider a strategy for getting the space from V((f.. By the online 

basis functions we mean basis functions that are computed during iterative process, contrary to offline basis 
functions that are computed before iterative process. The online basis functions are computed based on 
some local residuals for the current multiscale solution 

Let = td((^ + span{(/5} be the new approximate space that constructed by adding online basis € V) 

on the i-th coarse neighborhood uji and G Kns"''^ be the corresponding GMsFEM solution. 

We define which satishes 

ar(c”+\,i;) = r(/,i;) +m(c”„i;). 
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For the error we have 


„ l'r”+i - r' 
“t I’-ms F 


,ra+l „"+l _ 

semi’^ms ^semi 


=/7 rr’^+^ ^ - r 

iJ “'TV'-'ms ’'^ms ‘^semi/ v^semi ’'^ms 

= G (') — (T f ') — TTlf ^ 

“"T’V'-'ms ’'^ms '^semi/ V' */’'^ms '^semi/ ''‘'V'^ms’*^ms '^semi/ 


'semi ’ ms 


.n+i 'i 
semi/ 


AT 


— r('c”+^ - c"+M — A^r/Y fPhc”^^ - c”+M + c”+^ - c”+M'l 

~ ' V^ms ^semi/ ~ Z-^ ' ^VA^V^semi ^ms ) ^ ‘-ms Semi/v' 
i 

' 2 M* M A^ v-^ 2 V‘-semi ‘-ms / ^‘-ms ‘-S' 


< 


N 


^semi/l l"^* 


Using (25) we obtain 


ar{c 


n+1 

ms 


_ „"+l „n+l 

‘-semi’ ‘^ms 




semi/ — 



(26) 


The solution 


satishes 


11 n+l,(m+l) 
M‘-ms 


-Ct:ill^<lk-c 


n+1 

semi 


2 

T’ 


VueU++i. 


Taking v = 


+ aif, we have 


II n+l,(m+l)_ n+l||2 ^ 
11 '^ms '^semi Hr — 


II n+1,(m) 

I +ms 




2 

T 


II n+1,(m) 
M'-ms 


„n+l 112 


+ 2 a a^^{c 


n+l,(m) 

ms 


— C 


n+1 
semi ’ 


+++a“*((^,+. 


The last two terms in above inequality measure the amount of the reduction in error when the new basis 
function ip is added to the space 

For a = —1 and ip £ Vi is the solution of 


ar{p,v) = r{v), \/v£Vi. 


Then for - c”±^= we have 


11 n+l,(m+l )||2 ^ II n+1, (m) m 2 
I+semi Mr — I+semi Mr 


< iip”+i/'"b |2 
— I +semi I It 


2a-pic:^^^’^^^-c:^Z‘P)+r^{‘p) 

2a--(c((,+’(-),+ + 2r(/,u) + 2m(cS,+),+ + r(+ < 


2 

^ ■ 


To enhance the convergence and efficiency of the online adaptive GMsFEM, we consider enrichment on non¬ 
overlapping coarse neighbothoods. Let I C {1,2,..., TV} be the index set of some non-overlapping coarse 
neighborhoods. We define = 1^ -I- span{</Ji, i G /} and obtain 


II n+l,(m+l )||2 ^ 
11‘-semi Mr — 


II n+l,(m) 
P • 

11 semi 




2 

■ 


iGl 


Finally, we combine this with (26) 


and obtain 


II n+l,{m+l )||2 ^ II n+l,(m )||2 
I+semi Mr — I+semi Mr 




‘^r v'-semi ’ ^semi / 



< 



E.^i\\n\\l \ 

{c+^)eImi) 


II "+l,(m)||2 
I+semi Mr' 


We will find online basis functions p £ Vi to maximize the local resudial r" for current time level. 
Moreover, the required p is the solution of 


ar{p,v) = ri{v), Vu G U, 


(27) 
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where Ri{v)is the local residual that defined using 

r^(v) = T(f,v) - 

and ||ri||^ = ||(/3||^ according to the Riez representation theorem. 

For solution in each time level, we iteratively enrich our offline space by residual based online basis 


function. These basis functions are calculated using Equation (271 with zero Dirichlet boundary conditions 
and the residual norm ||r”||* provides a measure on the amount of reduction in energy error. 

For the construction of the adaptive online basis functions, we first choose 0 < 0 < 1, for each coarse 


neighborhood Wj, find the online basis tpi G Vi using equation (27). After compute the norm of local residuals 
and calculate rji 

Vf ■■= llr.llS, 

where ||ri||* = ||(/>i||r, then arrange them in descending order, i.e. r]1 > rj^ > ... > ? 7 ^. Then, choose the 
smallest k such that 

N k 

i=l i=l 

This implies that, for the coarse neighborhood u}j[j = 1, ...A:), we add the corresponding online basis to 
the original space 


6.1 Numerical results 

Next, we present numerical results for residual based online basis functions. We consider a similar problem 
as in the previous section with constant matrix-fracture properties in (§ and iteratively enrich the offline 
space by online residual basis functions in some selected time steps. Our coarse and fine grid setups are the 
same as in Section [4.1| Because re-generation of the matrix R is needed, when we add online basis function, 
we add them for some selected time steps. We note that, when we add new online basis functions, which 
are based on current residuals, we remove previously calculated online basis function and keep them till 
we update the online basis functions. It will save computational time if we have small size of coarse scale 
problem. 

In Table we present and errors. We consider three different cases. In the first case (we call 
it Case 1), online basis functions are added at the first time step and after that in every 30-th time step. 
In the second case (we call it Case 2), online basis functions are replaced at the first five consecutive time 
steps, and after that, the online basis functions are updated in every 30-th time step. In the third case (we 
call it Case 3), online basis functions are replaced at the first ten consecutive time steps, and after that, 
the online basis functions are updated in every 30-th time step. More updates initially helps to reduce the 
error due to the initial condition. As we mentioned that the offline space is important for the convergence, 
and we present the results for different number of initial offline basis functions per coarse neighborhood. We 
use multiscale basis functions from offline space as a initial basis functions. In Tablewe show errors when 
online basis functions are replaced at the first five consecutive time steps (as in Case 2), and afterwards, 
online basis functions are updated at lO-th, 20-th and 30-th time step. For our calculations, we use t^ax = 5 
years with r = 10 days. Calculations are performed in the coarse grid with 121 nodes for the case with 
constant matrix-fracture properties. We observe from this table the following facts. 

• Choosing 4 initial offline basis functions improves the convergence substantially. This indicates that 
the choice of the initial offline space is important. 

• Adding online basis functions less frequently (such as at every 30th time step) provides an accurate 
approximation of the solution. This indicates that the online basis functions can be added only at 
some selected time steps. 

Next, we would like to show that one can use online basis functions adaptively and use the adaptivity 
criteria discussed above. In Table we present results for residual based online basis functions with adap¬ 
tivity with 9 — 0.7. In Figure [T^ we show errors by time. We observe that applying adaptive algorithm can 
much reduce errors. 


20 





Figure 12: Fine scale solution (right), coarse-scale using 2 offline basis functions (middle) and coarse-scale 
after two online iteration for some time levels (left) for t=l, 3 and 5 year (from top to bottom) (constant 
matrix-fracture coefficients). For fine-scale solution size of problem is 6659. For 2 offline basis functions is 
242 and after two online iteration is 484 


7 Conclusions 

In this paper, we present a multiscale approach for shale transport in fractured media. Our approach uses an 
upscaled model in the form of nonlinear parabolic equations to represent the matrix that consists of organic 
and inorganic matter. The nonlinearities in the equation are due to the interaction of organic and inorganic 
matter. The interaction of nonlinear matrix and the fracture is represented by multiscale basis functions. We 
follow Generalized Multiscale Finite Element Method to extract the leading order terms that represent the 
matrix and the fracture interaction. Multiscale basis functions are constructed locally in each coarse region 
and they represent the interaction between the upscaled matrix and the fracture network. We show that 
our proposed approach can effectively capture the small-scale effects and the overall system can be modeled 
using a fewer degrees of freedom. Numerical results are presented. In some cases and some regions, the 
offline procedure is insufficient to give accurate representations of the solution, due to the fact that offline 
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DOF 
(# iter) 


Hi 

Moff = 1 

121 

17.136 

68.989 

242 (1) 

13.047 

43.662 

363 (2) 

7.275 

12.262 

II 

242 

3.975 

36.337 

363 (1) 

1.653 

6.705 

484 (2) 

0.889 

0.972 

II 

484 

0.651 

3.595 

605 (1) 

0.208 

0.307 

726 (2) 

0.171 

0.056 


DOF 
(# iter) 


Hi 

Moff = 1 

121 

17.136 

68.989 

242 (1) 

13.209 

42.777 

363 (2) 

6.603 

12.689 

II 

242 

3.975 

36.337 

363 (1) 

1.716 

7.841 

484 (2) 

0.546 

0.914 

II 

484 

0.651 

3.595 

605 (1) 

0.165 

0.313 

726 (2) 

0.105 

0.057 


DOF 
(# iter) 


Hi 

Moff = 1 

121 

17.136 

68.989 

242 (1) 

13.002 

42.091 

363 (2) 

6.125 

13.186 

II 

242 

3.975 

36.337 

363 (1) 

1.692 

8.709 

484 (2) 

0.449 

0.862 

II 

484 

0.651 

3.595 

605 (1) 

0.144 

0.318 

726 (2) 

0.076 

0.056 


Table 6: Convergence history using one, two and four offline basis functions {Moff = 1,2 and 4). We add 
online basis functions for every 30 time step and for iV—th first steps (Cases 1, 2, and 3). Left: = 1. 

Middle: N = 5. Right: N = 10. Here DOF for the last time step. 


DOF 
(# iter) 


Hi 

Moff = 1 

121 

17.136 

68.989 

242 (1) 

12.511 

41.477 

363 (2) 

5.910 

13.153 

Moff = 2 

242 

3.975 

36.337 

363 (1) 

1.624 

8.225 

484 (2) 

0.378 

0.934 

II 

484 

0.651 

3.595 

605 (1) 

0.126 

0.303 

726 (2) 

0.048 

0.033 


DOF 
(# iter) 


Hi 

Moff = 1 

121 

17.136 

68.989 

242 (1) 

12.925 

42.778 

363 (2) 

6.275 

12.705 

Moff = 2 

242 

3.975 

36.337 

363 (1) 

1.669 

8.034 

484 (2) 

0.474 

0.880 

II 

484 

0.651 

3.595 

605 (1) 

0.147 

0.305 

726 (2) 

0.080 

0.041 


DOF 
(# iter) 


Hi 

Moff = 1 

121 

17.136 

68.989 

242 (1) 

13.209 

42.777 

363 (2) 

6.603 

12.689 

II 

242 

3.975 

36.337 

363 (1) 

1.716 

7.841 

484 (2) 

0.546 

0.914 

II 

484 

0.651 

3.595 

605 (1) 

0.165 

0.313 

726 (2) 

0.105 

0.057 


Table 7: Convergence history using one, two and four offline basis functions {Moff = 1,2 and 4). We add 
online basis functions for every TVth time step and for first 5 steps Left: N = 10. Middle: N = 20. Right: 
N = 30. Here DOF for last time step 


computations are typically performed locally and global information is missing in these offline information. 
These phenomena occur locally and in some of these regions that are identified using the proposed error 
indicators, we need to develop online basis functions [8j. We discuss online basis functions and show that 
this procedure converges fast. 
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DOF 
(# iter) 

E iter 


Hi 

Moff = 1 

121 


17.136 

68.989 

242 (1) 

11 

13.209 

42.777 

363 (2) 

22 

6.603 

12.689 

Moff = 2 

242 


3.975 

36.337 

363 (1) 

11 

1.716 

7.841 

484 (2) 

22 

0.546 

0.914 

Moff = 4 

484 


0.651 

3.595 

605 (1) 

11 

0.165 

0.313 

726 (2) 

22 

0.105 

0.057 


DOF 

E iter 


Hi 

Moff = 1 

121 


17.136 

68.989 

243 

44 

4.082 

6.418 

381 

78 

2.681 

2.871 

0^ 

II 

to 

242 


3.975 

36.337 

376 

44 

0.441 

0.720 

504 

77 

0.376 

0.325 

II 

484 


0.651 

3.595 

635 

44 

0.110 

0.044 

737 

68 

0.098 

0.039 


Table 8: Convergence history using one, two and four offline basis functions {Moff = 1,2 and 4). We add 
online basis functions for every 30 time step and for first 5 steps Left: without space adaptivity. Right: with 
space adaptivity. Here DOF for last time step 
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